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Abstract 

A generalized flow solver using an implicit Lower-Upper(LU) diagonal decomposition based numerical 
technique has been coupled with three low-Reynolds number k-e models for analysis of problems with engi- 
neering applications. The feasibility of using the LU technique to obtain efficient solutions to supersonic 
problems using the k-e model has been demonstrated. The flow solver is then used to explore limitations and 
convergence characteristics of several popular two equation turbulence models. Several changes to the LU 
solver have been made to improve the efficiency of turbulent flow predictions. In general, the low-Reynolds 
number k-e models are easier to implement than the models with wall-functions, but require much finer near- 
wall grid to accurately resolve the physics. The three k-e models use different approaches to characterize the 
near wall regions of the flow. Therefore, the limitations imposed by the near-wall characteristics have been 
carefully resolved. The convergence characteristics of a particular model using a given numerical technique 
are also an important, but most often overlooked, aspect of turbulence model predictions. It is found that some 
convergence characteristics could be sacrificed for more accurate near-wall prediction. However, even this 
gain in accuracy is not sufficient to model the effects of an external pressure gradient imposed by a shock- 
wave/boundary-layer interaction. Additional work on turbulence models, especially for compressibility, is 
required since the solutions obtained with base line turbulence are in only reasonable agreement with the 
experimental data for the viscous interaction problems. 


Nomenclature 



A,B 

Convective Jacobians 

M 

Mach Number, (=u/(yRT) 1/2 ) 

a,b,g,h 

Constants in Van-Driest Transformation, 
Eqn 19, Eqn 20 

N 

Conservative/Non-conservative Transfor- 
mation Vector. 

D 

Difference Operator 

Pk 

Turbulence Production 

^el^e2 

Turbulence Model Constants, Table 1 

Pr 

Prandtl Number 


Eddy Viscosity Constant, (=0.09) 

P.P 

Pressure 

c f 

Skin Friction Coefficient 

Rt 

Turbulent Reynolds Number(=k 2 /ve) 

C, 

A Parameter in Eqn. 20 

Rk 

Near Wall Turbulence Reynolds Number 

E.F.Q 

Conservative Vectors 


(=k 1/2 y n /v) 

f H 

Wall Damping Function, Eqns 4, 5, 6 

R 

Gas Constant 

H 

Source Vector 

r 

Recovery factor(=Pr l/3 ) 

H 

Enthalpy 

T 

Temperature 

h 

Cell Volume 

T 

‘u 

Turbulence Intensity 

I 

Identity Vector 

u,v 

Mean Velocity Constants 

k 

Turbulent Kinetic Energy 

u T 

Friction Velocity (=(x w /p) 1/2 > 

( 

Length Scale, Eqn 22 

x,y 

Coordinates 

M n> M^ 

Implicit Viscous Terms 

y + 

Non-Dimensional Distance 


Symbols 

e Dissipation 

8 Boundary Layer Thickness 

y Maximum spectral radii (=MAX(X)) 

y, Ideal gas constant(= 1 .4) 

X Local spectral radii 

A Low Reynolds number terms 

|i Viscosity 

Eddy Viscosity(=C^k 2 /£) 
v Kinematic Viscosity 

r|,£ Body fitted coordinate system 

n Wake Constant(=0.50) 

o Viscous diffusion 

o Viscous Diffusion Coefficient, Table 1 

x Shear Stress 

Subscripts and superscripts 

a,b With respect to the Directions 

be Boundary Condition 

k With respect to Turbulent Kinetic Energy 

n pseudo time index 

t With respect to Turbulence 

v With respect to Viscous Terms 

A Transformed Variable 

oo Freestream 

0 At the inlet 

During the past several years, significant 
progress has been made in the development of new turbu- 
lence models and in the understanding of the turbulence 
model limitations. A majority of modelling efforts have 
been directed towards three classes of turbulence models 
with various degrees of physical complexity. It was previ- 
ously shown that the simplest turbulence models can be 
used to perform reasonably accurate combustor flowfield 
predictions. 1 This work also showed that a more accurate 
turbulence model may be needed to resolve the desired 
performance characteristics of the combustors. Therefore, 
prior to studying a very complex problem in a combustor 
flowfield, the turbulence model and the numerical tech- 
niques in non-reacting situations were studied and are 
described herein. 

The simplest of the turbulence models usually assume 
length scale behaviors and are known as the zero equation 
algebraic eddy viscosity models. These models utilize local 
equilibrium turbulent scale distributions and the mean 
shear information to complete the Boussinseq hypotheses 
necessary to close the averaged Navier-Stokes equations. 
The models in this class are widely used in Computational 


Fluid Dynamic(CFD) computations because of their sim- 
plicity. However, the major Achilles heel of these models 
is the limitation imposed by the equilibrium assumption 
used in the development of the length scale profiles. In 
flowfields involving non-equilibrium phenomena, such as 
separation, sudden removal or addition of pressure gradi- 
ent, injection, suction, wall temperature gradient and 
roughness, the predictions made using these models tend to 
be poor because the required model physics often exceed 
the equilibrium assumption. 2 

The second class of model retains the eddy viscosity/ 
Boussinseq approximations. However, these models use 
transport equations to find the required turbulent length and 
velocity scales of the flowfield to complete the eddy vis- 
cosity formulation and therefore these models are not 
limited to equilibrium situations. Unfortunately, a majority 
of the modeling parameters and the near-wall formulations 
are optimized only for simple flow situations. Therefore, 
the limitations of these models must be carefully explored 
in complex flow situations. There are many viable sugges- 
tions as to which quantities should be used to characterize 
the flowfield. The turbulent kinetic energy is the first 
dependent variable used in many of these models. How- 
ever, the choice of the second dependent variable required 
to determine the length scales varies. Some of the more 
successful choices are dissipation rate(e), 3 specific dissipa- 
tion rate(co), 4 and dissipative time scale(x). 

For the third class of turbulence model, the Reynolds 
stresses(Tjj) are used directly to close the mass(Favre)-aver- 
aged equations. This is usually called the second order 
closure model or the Reynolds stress model. The second 
order closure model formulation eliminates the isotropic 
assumption in which the eddy viscosity is aligned only with 
the mean shear. However, this gain in generality does not 
always translate directly into improved solution accuracy. 6 
Recently, several investigators have reported that complex 
flow situations can indeed be resolved using the Reynolds 
Stress model and very promising predictions have been 
reported by Launder-Shima 7 and Lee et al. 8 However, 
these works have also shown that the solution to full Rey- 
nolds-stress models can be achieved only through 
considerable computational effort which is undesirable in 
engineering analysis. 

The most wisely used and extensively developed of the 
two equation eddy viscosity models is the k-e model. In this 
formulation the length scale and the velocity scale needed 
to determine the eddy viscosity are computed from the 
modeled turbulent kinetic energy and dissipation rate trans- 
port equations. The k-£ model has been used in a wide 
variety of low Mach number problems. There are, however, 
several well known k-£ model deficiencies. A major defi- 
ciency of the model is its inability to handle no-slip 
surfaces. There are two techniques commonly used to 
resolve this deficiency. Extensive review of this deficiency 
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and analysis of the various model correction formulations 
have been conducted by Rodi 9 , Coakley and Huang, 10 and 
subsequently by Michelassi and Shih. 1 1 The first method 
replaces the turbulence transport equations in the near-wall 
regions with a known solution, usually based on the log- 
law assumption. This is usually known as the high-Rey- 
nolds number/wall-function approach. The second method 
is to modify the transport equations so that the singular 
behavior in the near wall region is removed and much sim- 
pler zero boundary conditions are satisfied. Only the low- 
Reynolds number models were studied within this work. 

A majority of the models previously discussed were 
developed for incompressible flows and have not been 
extensively tested for high Mach number flows. The com- 
pressible versions of these models, applicable in high Mach 
number flows, are just being introduced and potential prob- 
lems are being realized. 10 The early compressible 
turbulence models have been shown to be inadequate due 
to the lack of dilatational effects in the original formula- 
tions. Therefore, extensive efforts 7,8 are being made to 
construct a physically optimal turbulence model for high 
Mach number compressible flows. These efforts have 
resulted in considerable advances being made in the devel- 
opment of models for various dilatation terms, such as 
dilatation pressure and dilatation dissipation. 12,13,14 How- 
ever, these new modeling concepts have not yet been 
extensively tested for combustor flowfields. Therefore, in 
this work we have coupled the latest two-dimensional 
implicit finite volume numerical technique with the latest 
compressible two equation turbulence models, including 
the effects of dilatational terms, to develop a generalized 
CFD solver that can be economically used to solve a vari- 
ety of complex internal flow problems of engineering 
interest. 

The convergence characteristic is an important, but 
most often overlooked, aspect of any analysis using turbu- 
lence models. 10 In many instances, the numerical 
techniques will introduce unwanted transients, usually 
associated with numerical error, that can destabilize the 
numerics. This problem becomes more acute when conver- 
gence acceleration techniques are used. In some situations, 
local changes can translate into irreversible and non-phys- 
ical turbulence characteristics. A typical example of this 
situation occurs in the computation of the ffeestream 
region, where the turbulence problem should reduce to a 
natural decay, but more often the computed turbulent 
energy decays more rapidly than expected because of the 
numerical dissipation. This unexpected dissipation rate can 
cause the turbulent kinetic energy to decay to non-physical 
values, causing numerical instabilities to occur and conse- 
quently leading to convergence failure. There are a number 
of different suggestions for removing this limitation. The 
first is to limit the turbulence production during the itera- 
tion process to physically realistic levels. A second method 


is to start the solution process from a previously converged 
solution. Neither of these methods is very palatable in engi- 
neering applications. Recently, Huang and Coakley 15 
demonstrated that neither of these arrangements is neces- 
sary if the numerics has certain Total Variation 
Diminishing(TVD) properties and if the impending oscilla- 
tions are not allowed to turn negative. However, it is 
difficult to maintain TVD-like properties for all of the 
existing numerical schemes. Therefore, we have attempted 
to incorporate some of the proposed numerical behavior 
into the solvers. 

The focus of this work is to describe the development 
of an efficient LU solver for the low-Reynolds number k-e 
models and to address some of the key issues regarding the 
turbulence model physics and their implications in high 
speed internal flow predictions. 


I, Generalized TXirbulence Model Formulations 

The governing equations of a typical two equation 
model can be written in the conservative form: 

™+& F - F - )+ fr G - a -' i = H 

(i) 

where 
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The form of a ky and c ey are similar to above. The pro- 
duction, shear stress(T xy ), and normal stress(T xx ,T yy ) terms 
are defined as follows. 
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where the turbulence modeling constants vary from model 
to model and are shown in Table 1. All of the quantities 
appearing above _are mass(Favre 16 ) averaged, i.e. u = 
P M /p, and k = pMf>. The density(p) is conventionally 
averaged. 

In general, low-Reynolds number k-e models use mod- 
ified transport equations, where a fictitious dissipation(e) 
rather than the real dissipation(e) is used. The fictitious dis- 
sipation is usually defined as £ = £ - £wall- There are 
typically three different types of low-Reynolds number k-e 
models. All three employ additional wall damping to sat- 
isfy any additional near wall physics. The first type uses the 
flowfield characteristics, such as curvature of the turbulent 
kinetic energy(e wa n = 2v(/k/dy 2 ) 2 ), to satisfy the near wall 
asymptotic behavior, while at same time avoiding the sin- 
gularity caused by the no-slip boundary condition. The 
original k-e model formulation proposed by Launder and 
Sharma is typical of this type of formulation. The second 
type of low-Reynolds number k-e model, proposed by 
Chien, 18 uses the wall normal distance(e wa n = 2vk/y n 2 ) to 
develop the necessary near wall characteristics. The third 
type of low-Reynolds number formulation makes no 
attempt to correct the transport equation; rather, finite near- 
wall characteristics are used as boundary conditions. The 
low Reynolds number model proposed by Shih 19 is typical 
of the models that avoid the zero boundary condition and 
therefore the singularity. In practice, the low-Reynolds 
number k-e models are easier to numerically implement 
than a high Reynolds number model. However, the low- 
Reynolds number models also require much higher near 
wall grid resolution. Typical modeling constants of these k- 
e models are shown in Table 1. 


Models 

^el 

c e2 

c. 



Chien 

1.35 

1.80 

0.09 

1.0 

1.30 

Launder 

Sharma 

1.45 

1.92 

0.09 

1.0 

1.30 

Shih 

1.45 

2.00 

0.09 

1.0 

1.30 


Table 1 : Turbulence Model Parameters 


La Low-Reynolds Number k-e model Formulation of 
Launder and Sharma(1976) 

As previously discussed, the k-e models employ 
additional damping functions to satisfy the near wall 
physics. All of the models try to satisfy y n 2 behavior of the 
turbulent kinetic energy using slightly different damping 
functions along with the e wal j formulation. The Launder- 
Sharma model employs a near wall damping formulation 
that satisfies the invariance principle by using only near 
wall scalar characteristics. Therefore, the damping 
functions are independent of the grid coordinate system. 
Thus, for complex wall geometries additional modeling 
assumptions are not needed and the near-wall 
characteristics are uniquely defined. Unfortunately, this 
model is more sensitive to numerical error than Chien’ s 
model. Similar difficulties have been observed by 
Gerolymos. 18 The reason for these difficulties lies in the 
sub-layer region of a typical boundary-layer computation 
where the initial transient causes a large and unphysical 
production of turbulent energies which sometimes cannot 
be readjusted to the field. This difficulty was controlled 
using a limiter on the dissipations in the work of 
Gerolymos. However, we have found that such a limiter 
is not needed if a low CFL number is used to reduce the 
initial transients. The low-Reynolds number correction 
terms for this model are defined as follows, 


-2|i (V Jk ) 2 



(4) 


A = expf- 


-3.5 


^ (50 + R,)‘ 


I.b Low Reynolds Number k-e model Formulation of 
Chien(1976) 

A low-Reynolds number k-e model with Van- 
Driest type near-wall damping function and algebraic near 
wall formulation was proposed by Chien. 18 This low- 
Reynolds number formulation is not Gallean-invariant and 
is not independent of the grid coordinate system used. 
Therefore, a major disadvantage of this formulation is that 
for complex geometries the near wall characteristics may 
not be uniquely defined; for some cases, there may be wide 
range of possibilities. However, for the problems 
investigated, this model seems to have the best overall 
convergence behavior and numerical stability of the three 
models explored. Furthermore, since the test problems 
investigated are bounded by a single no-slip surface, lack 
of invariance was not a problem. The low-Reynolds 
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number correction terms for the Chien model are defined as 
follows: 


A = 



-2| x^exp (~0.5y + n ) 

y n 


(5) 

u = 1.0- exp (0.01 15)0 

where y n is the nearest surface normal distance. 

I.c Low-Reynolds Number k-e model Formulation of Shih 
and Lumley(1990) 

This model formulation involves a hybrid 
approach in which the near wall characteristics are 
modeled as boundary conditions and not as modified terms. 
The typical form of the low-Reynolds number damping 
terms are defined as follows: 


0.0 



( 6 ) 


/„ = VO +exp(a l R k + a 3 R k + a s R k )) 

where the damping constants aj, a 3 , and a 5 are defined as 
-1.5xl0" 4 , -l.OxlO" 9 and -5.0xl0" 10 , respectively. Note that 
the dissipation term in the turbulent kinetic energy equation 
has not been modified. However, this formulation uses the 
surface normal distance in the damping function, as does 
Chien’ s formulation. Therefore, this model does not satisfy 
the invariant principal; further, it is sometimes difficult to 
determine the proper surface normal distance in complex 
flow situations. However, because the normal distance is 
used only in an exponential term of the damping function 
formulation, this model is far less sensitive to the choice of 
the surface normal distance than the Chien formulation, 
where the surface normal distance is used directly as a 
dissipation correction term. 

I.d Compressibility Effects 

The effect of compressibility is an important 
consideration in developing a high speed turbulence model. 
Essentially, disregarding the dilatation effects will lead to 
under representation of the additional dissipation caused by 
volumetric or density fluctuations. Therefore, the effects of 


both the dilatation dissipation and pressure transport terms 
are usually modeled using the solenoidal dissipation and 
the turbulent Mach number(M t ) or density variance. 14 The 
simplest form of the proven modeling concept proposed by 
Zeman 12 and Sarkar 13 has been used to close the dilatation 
dissipation and the pressure dilatation effects. Since near- 
wall turbulence is considered in this study, a more 
appropriate form with correct near-wall limit suggested by 
Wilcox 21 has been implemented. 

z d = a (Mj - 0.06) eH (M t - 0.25) 


where H(>1.0) = 1.0 and H(<0) = 0.0. 


(7) 


II. Numerical Techniques 

Implicit LU finite volume numerical techniques 
were used to solve the Favre averaged Navier-Stokes 
equations along with the transport equations of turbulent 
kinetic energy and dissipation. This LU numerical 
technique was originally proposed by Yoon and 
Jameson. 22 The details of the Successive Symmetric Over 
Relexation(SSOR) scheme are well documented. 
Therefore, only the specifics dealing with the numerics of 
the turbulent transport equations are discussed here. One 
advantage of the LU diagonal decomposition technique is 
that the large CFL number marching ability of an implicit 
numerical schemes to be retained without the large block 
matrix operations. Furthermore, when the LU decomposed 
implicit operators are combined with a directionally 
sweeping procedure, the solution requires only a series of 
simple back substitutions. 

The approximate splitting method proposed by 
Jameson 22 based on the maximum spectral radii has been 
used in this work. This leads to Lower and Upper matrices 
which then can be solved using scalar inversions. A more 
rigorous LU formulation Steger and Warming(SW) flux- 
vector splitting procedure was also studied. In general, the 
splitting process reduces the complexity of the numerical 
operations in return for reduced convergence 
behavior(Hirsh 24 ). In order to regain some of the lost 
convergence speed and accuracy, the viscous terms have 
been added implicitly and the explicit side of the numerics 
was refined using an upwind technique. These two 
improvements further enhance the skin friction 
convergence. Typical convergence characteristics of these 
numerical techniques are compared in Figure 1 . The time 
history of the skin friction coefficient at the last wall 
surface station of a boundary-layer is shown in the top 
figure, and the total residual of the mean terms is shown in 
the lower figure. This figure shows that a reasonable 
convergence behavior can be achieved using the numerical 
techniques investigated. This figure also shows that the S W 
left-hand side and upwind differencing(UP) improve the 
convergence speed. Typically, approximately 500 to 1000 
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iterations based on the third significant figure of the skin 
friction predicted separate the three options shown in this 
figure. This figure shows that the SSOR technique with 
central differencing(CD) has the slowest convergence rate. 


Convergence History 



0.0 1000.0 2000.0 3000.0 

Figure 1: Typical Convergence Characteristics of the LU 
Numerics 


The turbulence transport equations are solved uncou- 
pled from the mean equations in order to reduce the 
numerical complexity. Therefore, the implicit operators 
needed to solve the additional two equations do not 
increase the rank of the overall implicit operator matrix for 
the mean solver. Thus, the time marching of the turbulence 
transport equations lags the mean equations marching by 
one iteration. This decoupling is not a major problem 
because the link between the mean turbulence equations is 
weak and occurs only through the diffusion terms. While 
the decoupling greatly simplifies the algorithm, some extra 
iterations may be required to resolve any sudden changes in 
the numerical solutions. Furthermore, it is possible that this 
can cause numerical instabilities during early stages of the 
computation, when rapid changes in the solution are 
expected. In order to minimize this effect, the source inter- 
action term has been implicitly formulated and the iteration 
process was started from CFL of 1.0 or less. 


III. Numerics 

The governing equations are transformed into the 
following form when the body fitted coordinate system is 
used: 


where 

Q = hQ H = hH 

E = h (£/ + Tl G) E v = h (£/ v + il x G,) 

G = h &F + Tl v G) Gy, - h ($F V + Tl G v ) 

( 8 ) 

The delta form of Equation 8 can be written in the 
following implicit for:, 

LS~ l UAQ = A tRHS 

L = (I + At(D\A + +D' ry B*-A-B-t)) 

U = (I+At(D\A +D^B-A + -B*)) 

S = (7 + Ar(A + +B-A-B )) 

(9) 

where the explicit right handside(RHS) is defined as the 
following, 

RHS = (D^(F - F v ) +D T1 (G-Gv)) +H 

and AQ is defined as Q n+1 -Q n . 

The source term(H) contains the effect of 
turbulent energy transfer between the dissipative scales and 
energy containing scales as well as some of the low 
Reynolds effects. Therefore, proper numerical treatment of 
these terms is essential in turbulent flow solver 
development. In near-wall flow situations, this term tends 
to get extremely stiff and therefore in most cases the 
implicit representation of this term is essential in 
maintaining numerical stability. This difficulty is further 
extenuated when the turbulence model is expected to 
handle a rapid exchange in turbulent energy between the 
dissipative scales and the turbulent energy scales due to 
externally generated gradients such as shock-waves. 
Similar observations have been made by other 
investigators. 15 ’ 20 However, the exact formulation of the 
source flux Jacobian term varies from study to study. 
Initially, we have used the following 2x2 matrix for the 
source Jacobian on the left hand side of the Equation 9: 


T-. = 
u 


r k e 

- ( — ) 


Wk 

pe 




Huang and Coakley 15 showed that the source 
Jacobian should be formulated so that the overall numerical 
scheme maintains a positive left hand side. We have found 
from our experience that this particular point is not 
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essential for the SSOR technique because of the diagonal 
dominance generated by the maximum spectral radii. 
However, for the S W technique, maintaining a positive left 
hand side is more important because the local spectral radii 
do not always dominate. In order to maintain a positive left 
handside, only the negative diagonal terms of T z y were used 
along with the SW algorithm. 

m.a LU-SSOR 

In the Symmetric Successive Over Relaxation 
numerical technique, the global eigenvalue, directional 
splitting, and LU decompositions are used to set-up the left 
hand side(LHS) of Equation 8 such that only scalar 
inversions are required to solve the equation. The flux 
Jacobians A, + A,~ B, + and B' are constructed such that the 
eigenvalues of A + and B + are positive and the eigen values 
of A“ and B" are negative. The convergence is achieved 

through a four step sweeping process: 

./ y A y B At ~+ w 

A Qij = Af((/+Af(^ + ^)))«ffS+^A Wf yAa.i 

y A y B ~ n 

(/+A,( 4 + ^ )_A ^ )Aa; = 

* / At .+ 4 ~ // 

A Qij+ 


K = 


u o 
0 u 


,± <wai> 

A = 


and as well as B + , B' are similarly defined. 


( 12 ) 


When a directional splitting process similar to the 
LU-SSOR is used. Equation 10 reduces to the following 
four step procedure: 

At . * + M k N 


A = AtRHS + ^ (A + + ^ AqIij 


A£ At| A^ 2 At) 


-Af7 \AQ" = A q'j 


At .♦ M n N 

+ -r -(S +-T 5 - 

At| Ar( 


) ® i,y - 1 A Q i,j - 1 




Ul . IbI 


i e» = l' +i, ^ + k^ll 4C v 


- // 



M k N 

"AT 


) 


^ /// 

A<2i + w 

i + i.y 


./// pii At 7 - pin 
AQij = A Qij — + \ jAQi + \ j 


Ul \b\ M t\ N W « W// 

/+ Af | ^ JJ AQ-,j = AQ^ 


A^ ATI a% AT) 2 


y A y B ~ ~ in At 

(/+A,( 4 + ^i ))A<2i2 = AQiJ ~A^ BiJ+lAQi - J+i 

(II) 

where the split flux Jacobians are defined as, 

A = A±y a / ) 

and B+, B- are similarly defined. Here, the spectral radii are 
corrected by the viscous terms M^N / A^ 2 and M^N/Aq 2 . 

m.b LU-SW 

In order to improve the convergence of skin 
friction prediction, flux vector splitting procedure of Steger 
and Warming 23 was also studied. This decomposition is 
achieved by using local eigenvalues based on local 
characteristic wave speeds. This approach yields simple 
rank two matrix operations for the de-coupled numerical 
formulation used in this study. The explicit terms on the 
right hand side of Equation 8 were differenced using the 
van Leer’s flux- vector upwind differencing technique to 
maintain numerical accuracy. Typical Jacobian matrices 
split according to the local wave speeds are as follows 


Af M^N 

—£r) A(2iy+i 
At) At] i+1 


(13) 

Here the viscous terms in thin-layer form have been added 
to the left handside to maintain stability in the near wall 
regions of the computation. 

III.c Van Leer’s MUSCL Upwinding Technique 

In order to obtain upwind numerical accuracy, the 
flux vector of the inviscid terms in the RHS are evaluated 
based on the Monotone Upwind-Centered Scheme for 
Conservative Law (MUSCL) 24 approach given by the 
following: 

^ = (F- (fi-) +r (< 2 + )) t - (F* (Q‘) +r (G + ))« 

F+ = ^pJ7RT(M±l) 2 Q 

(14) 

where non-conservative variables Q(k,e) are evaluated 
using the MINMOD flux limiter approach proposed by van 
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Leer(Hirsh 22 ). We have used a second order scheme for 
all of the computations. 

IV. Boundary and Initial Conditions 

Zero gradient conditions were assumed at the 
freestream boundaries of the computation. At the wall 
surfaces, zero boundary conditions were used in the low- 
Reynolds number k-e model of Chien 18 and Launder- 
Sharma. 17 Boundary conditions based on the near-wall 
Kolmogorov scales were used in the k-e model proposed 
by Shih. 19 The wall turbulent kinetic energy and 
dissipation are defined as follows: 

For Chien and Launder-Sharma Models, 

kbc = 0.0, and e^ = 0.0 

For Shih and Lumley Model, 

pkbc =0.25p w U t 2 , and pe^ = 0.25 1 p w 2 U x 4 /|i w 


( 17 ) 

This generalized transformation reduces to a 
simpler form derived by White 26 when flow is assumed to 
be adiabatic(T aw /T w = 1.0). Furthermore, the velocity 
profile is determined from this relationship given only the 
initial skin friction coefficient, boundary-layer thickness, 
and the temperature ratio. The thermal profile is developed 
from the modified Crocco-Cohen relationship, which is 
general enough to satisfy both cold wall and adiabatic wall 
boundary conditions. The Crocco-Cohen relationship of 
total enthalpy ratio has been shown by Wallace to be: 


(H-HJ 

: n e -Hj 



The inlet velocity profile was developed from 
the initial value of the skin friction coefficient, 
boundary-layer thickness and the law of the wall 
assumption. The Van-Driest transformation 25 was used 
to determine the velocity profile from the assumed log- 
law profile. The thermal profile was developed using a 
typical Crocco’s parabolic relationship. The turbulence 
kinetic energy and dissipation profiles were developed 
from the velocity profile, along with the equilibrium 
constraints and the assumed length scale profile: 

Typical law of the wall formulation can be 
written in the following non-dimensional form. 


y *-o3T" ,<y * >+5, + § H ' 


~y~ 

_5_ 


(15) 

where the function W(y/8) is the wake correction. The 
transformation proposed by Van Driest is then used to 
find the actual velocity profile which includes the effect 
of density from the log-law. A typical transformed 
equation for the arbitrary temperature ratio can be 
written in the following form: 



+ 2a 
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b = + a- 1 
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(18) 

where the parameter(C t ) ranges from 0, for Crocco- 
Busemann linear formulation, to 1-H e /H w , for a parabolic 
formulation which corrects for heat dissipation(Pr=1.0). 
This relationship can be further simplified by expanding 
the total enthalpy expression to the static temperature 
ratios. This simplification leads to the following expression 
for the non dimensional static temperature. 



g = (1 -c,) 
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Thus, once the velocity profile is developed from 
the previously derived velocity formulation the thermal 
profile can be completely developed from the above 
expression. Finally, the turbulent kinetic energy and 
dissipation rate profiles were developed from the 
equilibrium(P k =pe) assumption and the freestream 
condition constraints. 


k = 





* 3/2 
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+ e 


y 
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where the length scale is assumed to be 


( 20 ) 


l = min (0.4y, 0.0858) 


( 21 ) 
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The freestream turbulent kinetic energy and 
dissipation are determined from the following simple 
relationships, which maintain reasonable freestream 
turbulence and eddy viscosity levels: 



( 22 ) 


V. Test Cases 

Several supersonic turbulent boundary-layers 
with or without pressure gradients were analyzed to 
validate numerical implementation of the baseline 
turbulence models. Figure 2 shows typical Mach 2.87 inlet 
velocity and temperature profiles. These profiles have been 
developed by the method described previously to satisfy 
the experimental Reynolds number and the boundary-layer 
thickness. Initially, the profiles were computed a location 
upstream of the measurement station. Then the zero- 
pressure gradient boundary-layer problem was solved to 
find fully developed boundary-layer velocity and 
temperature profiles. These profiles are then used as the 
inlet conditions for the Navier-Stokes/turbulence-model 
computations. Typically, a boundary thickness of 0.025m 
and skin friction coefficient of 0.0013 have been used to 
develop the inlet profiles. A turbulence intensity of 1 
percent and an adiabatic wall have been assumed. 



A Mach 2.87 flow with freestream Reynolds 
number of 6.3xl0 7 /m was used to match the test condition 
of experimental studies of Smits et al. 28 ’ 29 A typical grid 
system used to resolve a zero-pressure gradient turbulent 
boundary layer is shown in Figure 3. The overall height of 
the grid system is usually three boundary-layer thicknesses. 
The length of the computational domain varied from 0.25 
m to 1 .00m. Uniform grid distribution is usually used in the 
streamwise direction. Typically, an 81 by 63 grid system 
was used for the shock-wave/boundary-layer interaction 
prediction and a 101 by 63 grid system was used in the flat 
plate boundary-layer predictions. The vertical spacing of 
this grid system has been exponentially stretched to 


maintain a near-wall spacing of approximately 1.0 in non- 
dimensional coordinate(y + ). Approximately, 45 of 63 
points are located in a boundary-layer region. The zero 
gradient boundary condition has been used in the exit and 
freestream plane of the grid system. The inlet profiles of 
velocity, temperature and pressure are fixed. 


Figure 3: A Typical Grid System for Zero Pressure 
Gradient Turbulent Boundary-layer Predictions 

Typical convergence behavior of the turbulent 
quantities along with the skin friction coefficient time 
history, obtained using a constant CFL variable time step 
marching technique, is shown in Figure 4. Convergence 
was determined by a two to three order reduction in the 
global residual of all quantities and reduction of the skin 
friction level to the third significant figure. This figure once 
again shows that the SSOR convergence is slower than the 
SW convergence. Although the SSOR numerical technique 
is stable to infinitely large CFL numbers, only values near 
one have been used to yield more meaningful comparisons 
to the S W flux vector splitting technique. This figure shows 
that SSOR technique with central difference RHS required 
more than 2000 iterations as compared to the 
approximately 1000 iterations required for the SW 
technique. 

io 4 

io J 
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1 ,0 " 
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Figure 4: Typical of Convergence Characteristics 



Figure 5 shows the typical skin friction coefficient 
comparisons between the model predictions and the 
experimental data reported by Smits et al. Although all of 
the model computations were started from the same inlet 
condition, all three models responded differently. The 
differences in the initial skin friction level can be as much 
as 60 percent. However, the skin friction differences relax 
into a more acceptable level down stream of the inlet, 
where the variation in the skin friction level reduced to 
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approximately 15 percent. The transition like behavior 
shown in this figure has been forced by not including the 
production contribution of the turbulence model in the first 
few stations of computation which force the solution to be 
laminar. Figure 5 also shows that Chien’s model prediction 
of the skin friction is the lowest of the three model 
predictions. 



Figure 5: Skin Friction Distribution Mach 2.87 Turbulent 
Boundary-Layer 

The law of wall profile predictions are compared 
in Figure 6. All profiles shown in Figure 6 have been 
transformed into equivalent velocity form using the Van- 
Driest transformation to eliminate the effects of variable 
density. Although the prediction made using the Shih- 
Lumley and Launder- Sharma models seem to be in closest 
agreement with the experimental data, all of the predictions 
are in reasonable agreement with both the experimental 
data and an earlier Reynolds Stress model prediction. 8 The 
deviation shown is mostly due to the discrepancy in the 
skin friction predicted by each of the models. 

Law of Wall Profile 



Figure 6: Typical Log-Law Profile Comparisons 


Typical profiles of the normal Reynolds 
Stresses(i xx x yy ) are shown in Figure 7. The Reynolds 
Stresses in this figure are normalized by the friction 
velocity. Note that the agreement between the experimental 


data and the prediction is only reasonable. Furthermore, 
because of the isotropic assumption, k-e model prediction 
of T xx is equal to T yy . This is not supported by the 
experimental data. This is consistent with the earlier 
observation of Lee. 9 Either an anisotropic eddy viscosity k- 
e model or a full Reynolds Stress model would be required 
to correctly model the normal stresses. Despite this 
difficulty, the other parameters of interest, such as skin 
friction coefficient, are predicted correctly. This is in part 
due to the damping factors employed by the k-e model 
which mask such deficiencies. 



Figure 7: Non-Dimensional Normal Stress(x xx T yy ) 


Figure 8 compares the measured kinematic eddy 
viscosity to the profile predicted by the k-e models. 
Although the data scatter is large, the predicted profiles are 
in reasonable agreement with the measured data. The major 
difference occurs in the outer region of the profile which 
seems to indicate that the freestream dissipation used may 
be too large. Furthermore, the peak value of eddy viscosity 
predicted and the high Mach number data are in reasonable 
agreement. 



0.0 1.0 
y/delta 


Figure 8: Typical Eddy Viscosity Profile 

A more complex problem of interest is the 
viscous-inviscid interaction caused by a shock-wave 
impinging on a turbulent boundary-layer. As shown later, 
the turbulence models seem unable to correctly handle 
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near-wall relaxation behavior caused by external pressure 
gradients. This is an important flow feature that contributes 
significantly to the flow losses in complex high speed 
combustors 30 and to other external aerodynamic 
parameters. 31 Therefore, accurate modeling of this 
problem is a primary requirement for any generalized CFD 
method developed for internal flow analysis. In order to 
study this class of problem in more detail and to perform a 
more rigorous test of the turbulence model-numerics 
combination developed, several cases of the viscous 
interaction problem experimentally studied by Smits et al. 
have been computed. In this flow configuration the shock- 
wave is generated by a wedge ramp The shock-wave then 
interacts with the incoming turbulent boundary-layer, 
leading to separation at the higher pressure ratio 
conditions. 

Turbulence model predictions of 8 and 24-degree 
ramp angle conditions are compared with the experimental 
data of Smits et al. 28,29 For the 8 degree case, the rise in 
pressure was high enough to cause only a thinning of the 
near wall boundary-layer, and not a separation. However, 
for the 24 degree case, the turbulent boundary-layer near- 
wall is fully separated. The separation makes the 24 degree 
cases a more difficult test for the turbulence models. A 
typical grid system used to model this flowfield is 
illustrated in Figure 9. The incoming boundary-layer for all 
of the predictions have been generated using Chien model. 
Typically, a flat plate region 0.05m and 0. 15m ahead of the 
comer has been resolved in the 8 degree case and 24 degree 
case, respectively. The wedge length is 0.2m, and a 
uniform grid distribution is used in the streamwise 
direction. The total height of the computational domain is 
three times the boundary-layer thickness. The vertical 
spacing of this grid system has been exponentially 
stretched to maintain a near-wall spacing of approximately 
1.0 in non-dimensional coordinate^*) at the inlet of the 
computational domain. 



Figure 9: Typical Flat plate-Wedge Grid System 


Figure 10 compares the predicted wall pressure 
and skin friction distributions and the experimental data of 
Smits et al. for the 8 degree case. This figure shows that the 
actual turning angle of the incoming turbulent boundary- 
layer is reasonably captured by all of the k-e models, which 
results in close agreement of wall pressure predicted and 
measured. However, the prediction of the skin friction is 
more difficult. The relaxation of the skin friction is over 


predicted by Chien model and is under predicted by the 
Launder-Sharma model. 



Figure 10: Surface Pressure and Skin friction Coefficient 


Distribution for 8 degree Ramp case 

Figure 1 1 compares the typical wall pressure and 
skin friction distributions predicted and the experimental 
data of Smits et al. for the 24 degree case. The presence of 
a strong shock-wave causes the boundary-layer to separate 
and to form a separation bubble near the comer of the 
compression ramp. The separation bubble then causes the 
inviscid turning well ahead of the wedge which leads to the 
two step wall pressure distribution. This early turning 
causes flow to developed a "lambda" shock structure with 
a pressure "plateau." 31 This flow structure has been 
predicted by all of the k-e models. The total rise in wall 
pressure has also been correctly predicted by all of the 
models. However, the predicted size of the bubble and the 



Figure 11: Surface Pressure and Skin friction Coefficient 


Distribution for 24 degree Ramp case 

strength of the wall skin-friction recovery vary among the 
models. Both the Launder-Sharma model and the Chien 
model over predicted the downstream skin-friction 
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recovery and the size of the separation bubble. The Shih- 
Lumley model prediction of the separation bubble is in best 
agreement with the experimental data. 

Figure 12 shows the effect of the dilatation 
dissipation on the 24-degree ramp flow prediction. Here, 
the k-£ model of Launder-Sharma is used as the base line 
model. The dilatational dissipation of Equation 7 causes 
the separation bubble to grow and wall skin-friction 
recovery to weaken(i.e. £=£+£ d ). Further addition of the 
turbulent kinetic energy to the thermal energy balance 
causes even further increase in the size of the separation 
bubble predicted(i.e. E^yT+U^+k). This observation is in 
agreement with Huang and Coakley. 15 Furthermore, 
Figure 12 shows that the compressibility terms are not the 
cause of the large separation bubble predicted by the 
baseline k-e model. Therefore, further study is needed to 
develop additional turbulence model corrections to 
improve separation bubble size predictions. 



X(m) 

Figure 12: Effects of Compressibility on the Prediction of 
24-degree Ramp case 


VI. Concluding Remarks 

A generalized flow solver using Lower-Upper 
diagonal decomposition techniques has been successfully 
coupled with three low-Reynolds number k-e models for 
engineering analysis. The limitations of the models and 
convergence characteristics in several supersonic flow 
situations have been explored. We have shown, that the 
combination of the LU based implicit numerical method 
and the low-Reynolds number k-e model can be used 
efficiently in studying problems ranging from simple flat 
plate boundary characteristics to complex shock-wave/ 
boundary-layer interaction. In this study we have also 
shown that an SSOR marching process is slowest in 
obtaining a steady state solution; and upwinding the 
inviscid terms tend to improve the efficiency. Furthermore, 
using the Steger- Warming flux vector formulation with the 
explicit terms upwinded seems to be optimal for choice for 
an LU based turbulence model solver. All three turbulence 


models are able to predict a zero pressure gradient 
supersonic turbulent boundary-layer with reasonable 
accuracy. However, the more complex problem of a shock- 
wave/boundary-layer interaction is handled less adequately 
by these models. The wall skin friction recovery behaviors 
predicted are drastically different among the models. It is 
found that for the 8 and 24 degree shock-wave/boundary- 
layer interaction studied, Launder-Sharma model and Shih- 
Lumley model predictions are closest to the experimental 
data. It was also found that the optimal CFL number which 
is allowed for both these model is approximately 1 .0 or less 
for the Launder-Sharma and Shih and Lumley model while 
the allowable CFL number for the Chien model can be up 
to 5. The effects of compressibility, in the form of 
dilatation dissipation, increase the separation bubble size 
and reduce the wall skin-friction level. 
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